function dif_B = PMSquFunc(x,y,z,xd,yd,zd,flag)
%integrand of flux of a cylindrical PM

R1 = sqrt((xd-x).^2+(yd-y).^2+(zd-z).^2);

switch (flag)
    case 1 
        dif_B = (xd-x)./R1.^3;
    case 2
        dif_B = (yd-y)./R1.^3;
    case 3
        dif_B = (zd-z)./R1.^3;
end
        
